!############################# Change Log ##################################
! 5.0.0
!
!###########################################################################
!  Copyright (C)  1990, 1995, 1999, 2000 - All Rights Reserved
!  Regional Atmospheric Modeling System - RAMS
!  Mission Research Corporation / *ASTeR Division
!###########################################################################

subroutine exevolve(mzp,mxp,myp,ia,iz,ja,jz,izu,jzv,mynum,edt,key)

use mem_tend
use mem_basic
use mem_grid
use mem_stilt
!use node_mod

implicit none

character(len=*)            :: key
integer                     :: mzp,mxp,myp,ia,iz,ja,jz,izu,jzv,mynum,i,j,k
real                        :: edt


if(key.eq.'ADV')then
! Initialization
call thvlastzero(mzp,mxp,myp,ia,iz,ja,jz,stilt_g(ngrid)%thvlast(1,1,1))

! Calculate advective term
  call exadvlf(mzp,mxp,myp,ia,iz,ja,jz,izu,jzv &
    ,grid_g(ngrid)%rtgu     (1,1)   ,grid_g(ngrid)%fmapui    (1,1)    &
    ,grid_g(ngrid)%rtgv     (1,1)   ,grid_g(ngrid)%fmapvi    (1,1)    &
    ,grid_g(ngrid)%f13t     (1,1)   ,grid_g(ngrid)%f23t      (1,1)    &
    ,grid_g(ngrid)%rtgt     (1,1)   ,grid_g(ngrid)%fmapt     (1,1)    &
    ,grid_g(ngrid)%dxt      (1,1)   ,grid_g(ngrid)%dyt       (1,1)    &
    ,basic_g(ngrid)%uc      (1,1,1) ,basic_g(ngrid)%dn0u     (1,1,1)  &
    ,basic_g(ngrid)%vc      (1,1,1) ,basic_g(ngrid)%dn0v     (1,1,1)  &
    ,basic_g(ngrid)%dn0     (1,1,1) ,basic_g(ngrid)%wc       (1,1,1)  &
    ,tend%pt                (1)     ,basic_g(ngrid)%pc       (1,1,1)  )
! Calculate compression term
  call excondiv(mzp,mxp,myp,ia,iz,ja,jz,izu,jzv                       &
      ,basic_g(ngrid)%uc    (1,1,1) ,basic_g(ngrid)%vc       (1,1,1)  &
      ,basic_g(ngrid)%wc    (1,1,1) ,basic_g(ngrid)%pc       (1,1,1)  &
      ,tend%pt       (1)            ,grid_g(ngrid)%dxt       (1,1)    &
      ,grid_g(ngrid)%dyt    (1,1)   ,grid_g(ngrid)%rtgt      (1,1)    &
      ,grid_g(ngrid)%rtgu   (1,1)   ,grid_g(ngrid)%rtgv      (1,1)    &
      ,grid_g(ngrid)%f13t   (1,1)   ,grid_g(ngrid)%f23t      (1,1)    &
      ,grid_g(ngrid)%fmapt  (1,1)   ,grid_g(ngrid)%fmapui    (1,1)    &
      ,grid_g(ngrid)%fmapvi (1,1)   )
! Put theta_v from last timestep into memory
  call fthvlast(mzp,mxp,myp,ia,iz,ja,jz,stilt_g(ngrid)%thvlast(1,1,1) &
      ,basic_g(ngrid)%theta (1,1,1) ,basic_g(ngrid)%rtp      (1,1,1)  &
      ,basic_g(ngrid)%rv    (1,1,1) )
elseif(key.eq.'THV')then
! Calculate heating term
  call exheat(mzp,mxp,myp,ia,iz,ja,jz,izu,jzv,mynum,edt,stilt_g(ngrid)%thvlast(1,1,1))
endif

return
end

![ML
subroutine thvlastzero(mzp,mxp,myp,ia,iz,ja,jz,thvlast)
  implicit none
  integer                                  :: mzp, mxp, myp, ia, iz, ja, jz
  real, dimension(mzp,mxp,myp)             :: thvlast
  thvlast(1:mzp,ia:iz,ja:jz)=0.0
  return
end subroutine thvlastzero
!ML]



!====================================================================
subroutine fthvlast(mzp,mxp,myp,ia,iz,ja,jz,thvlast,theta,rtp,rv)

use micphys, only: level


implicit none
integer                     :: mzp, mxp, myp, k, i, j, ia, iz, ja, jz
real,dimension(mzp,mxp,myp) :: thvlast, theta, rtp, rv

if(level > 0)then
  do j=ja,jz
    do i=ia,iz
      do k=1,mzp
        thvlast(k,i,j)=theta(k,i,j)*(1.0+1.61*rv(k,i,j))/(1.0+rtp(k,i,j))
      enddo
    enddo
  enddo
else
  do j=ja,jz
    do i=ia,iz
      do k=1,mzp
        thvlast(k,i,j)=theta(k,i,j)
      enddo
    enddo
  enddo
endif

return
end

!=========================================================================
subroutine exheat(mzp,mxp,myp,ia,iz,ja,jz,izu,jzv,mynum,edt,thvlast)

!use node_mod
use mem_basic
use mem_tend
use mem_grid

implicit none

integer :: mzp,mxp,myp,ia,iz,ja,jz,izu,jzv,mxyzp,mynum,i,j,k

real :: edt

real, dimension(mzp,mxp,myp) :: thvtend,thvadv,thvlast

!  Calculate theta_v tendency
call thetvtend(mzp,mxp,myp,ia,iz,ja,jz,izu,jzv                      &
   ,basic_g(ngrid)%theta    (1,1,1)  ,basic_g(ngrid)%rtp    (1,1,1) &
   ,basic_g(ngrid)%rv       (1,1,1)  ,thvtend,edt,thvlast)

!  Calculate advective of theta_v
call exthvadv(mzp,mxp,myp,ia,iz,ja,jz,izu,jzv,mynum,thvadv,edt)

!  Basically sum the above two contributions
call exhtend(mzp,mxp,myp,ia,iz,ja,jz                                 &
   ,basic_g(ngrid)%pi0      (1,1,1)  ,basic_g(ngrid)%pc      (1,1,1) &
   ,basic_g(ngrid)%rtp      (1,1,1)  ,basic_g(ngrid)%theta   (1,1,1) &
   ,basic_g(ngrid)%rv       (1,1,1)  ,tend%pt                (1)     &
   ,thvtend,thvadv)

return
end


!===========================================================================

subroutine thetvtend(mzp,mxp,myp,ia,iz,ja,jz,izu,jzv,theta,rtp, &
  rv,thvtend,edt,thvlast)

use micphys, only: level
use mem_grid, only: time

implicit none

integer :: mzp,mxp,myp,ia,iz,ja,jz,izu,jzv,mxyzp,i,j,k

real :: edt,edti

real, dimension(mzp,mxp,myp) :: theta,rtp,rv,thvtend
real, dimension(mzp,mxp,myp) :: thvlast

if(time > 0.1)then
    do j=ja,jz
      do i=ia,iz
        do k=1,mzp
          thvtend(k,i,j)=0.0
        enddo
      enddo
    enddo
  return
endif

edti=1.0/edt

if(level.gt.0)then
  do j=ja,jz
    do i=ia,iz
      do k=2,mzp-1
        thvtend(k,i,j) = (  theta(k,i,j) * (1.0+1.61*rv(k,i,j))  &
          / (1.0 + rtp(k,i,j))   -  thvlast(k,i,j)  )  * edti
      enddo
    enddo
  enddo
else
  do j=ja,jz
    do i=ia,iz
      do k=1,mzp
        thvtend(k,i,j) = ( theta(k,i,j) - thvlast(k,i,j) ) * edti
      enddo
    enddo
  enddo
endif  

return
end

!=======================================================================

subroutine exthvadv(mzp,mxp,myp,ia,iz,ja,jz,izu,jzv,mynum,thvadv,edt)

use mem_basic
use mem_scratch
use mem_grid
!use node_mod

implicit none

integer :: mzp,mxp,myp,ia,iz,ja,jz,izu,jzv,i,mxyzp,mynum,j,k

real :: dtlto2,edt

real, dimension(mzp,mxp,myp) :: thvadv,thetav,srthtv

mxyzp=mxp*myp*mzp

dtlto2 = .5 * edt

call prep_vt3d(mxyzp,dtlto2 &
   ,scratch%vt3da         (1)          , scratch%vt3db      (1)     &
   ,scratch%vt3dc         (1)          , basic_g(ngrid)%uc  (1,1,1) &
   ,basic_g(ngrid)%vc     (1,1,1)      , basic_g(ngrid)%wc  (1,1,1) )

call prep_thetv(mzp,mxp,myp,ia,iz,ja,jz &
   ,basic_g(ngrid)%theta  (1,1,1)      , basic_g(ngrid)%rtp(1,1,1) &
   ,basic_g(ngrid)%rv     (1,1,1)      , thetav)

call fa_preptc(mzp,mxp,myp &
   ,scratch%vt3da  (1)     ,scratch%vt3db  (1)     ,scratch%vt3dc  (1)      &
   ,scratch%vt3dd  (1)     ,scratch%vt3de  (1)     ,scratch%vt3df  (1)      &
   ,scratch%vt3dh  (1)     ,scratch%vt3di  (1)     ,scratch%vt3dj  (1)      &
   ,scratch%vt3dk  (1)     &
   ,basic_g(ngrid)%dn0    (1,1,1) ,basic_g(ngrid)%dn0u   (1,1,1)  &
   ,basic_g(ngrid)%dn0v   (1,1,1) ,grid_g(ngrid)%rtgt    (1,1)    &
   ,grid_g(ngrid)%rtgu    (1,1)   ,grid_g(ngrid)%rtgv    (1,1)    &
   ,grid_g(ngrid)%fmapt   (1,1)   ,grid_g(ngrid)%fmapui  (1,1)    &
   ,grid_g(ngrid)%fmapvi  (1,1)   ,grid_g(ngrid)%f13t    (1,1)    &
   ,grid_g(ngrid)%f23t    (1,1)   ,grid_g(ngrid)%dxu     (1,1)    &
   ,grid_g(ngrid)%dyv     (1,1)   ,grid_g(ngrid)%dxt     (1,1)    &
   ,grid_g(ngrid)%dyt     (1,1)   ,mynum)

do i=1,mxp
  do j=1,myp
    do k=1,mzp
      srthtv(k,i,j)=thetav(k,i,j)
    enddo
  enddo
enddo

call fa_xc(mzp,mxp,myp,ia,iz,1,myp,thetav,srthtv &
   ,scratch%vt3da  (1)     ,scratch%vt3dd  (1)     ,scratch%vt3dg  (1)      &
   ,scratch%vt3dh  (1)     ,scratch%vt3di  (1)     ,mynum)

if (jdim == 1)  &
  call fa_yc(mzp,mxp,myp,ia,iz,ja,jz,thetav,srthtv  &
      ,scratch%vt3db  (1)     ,scratch%vt3de  (1)  ,scratch%vt3dg  (1)      &
      ,scratch%vt3dj  (1)     ,scratch%vt3di  (1)  ,jdim,mynum)

  call fa_zc(mzp,mxp,myp,ia,iz,ja,jz,thetav,srthtv  &
      ,scratch%vt3dc  (1)     ,scratch%vt3df  (1)  ,scratch%vt3dg  (1)      &
      ,scratch%vt3dk  (1)     ,vctr1,vctr2,mynum)

do j=ja,jz
  do i=ia,iz
    do k=2,mzp-1
      thvadv(k,i,j)=0.0
    enddo
  enddo
enddo

call advtndc(mzp,mxp,myp,ia,iz,ja,jz,thetav,srthtv,thvadv,edt,mynum)

do j=ja,jz
  do i=ia,iz
    do k=2,mzp-1
      thvadv(k,i,j)=-1.0*thvadv(k,i,j)
    enddo
  enddo
enddo

return
end


!==========================================================================
subroutine prep_thetv(mzp,mxp,myp,ia,iz,ja,jz,theta,rtp,rv,thetav)
use micphys, only: level
implicit none

integer :: mzp,mxp,myp,k,i,j,ia,iz,ja,jz
real, dimension(mzp,mxp,myp) :: theta,rtp,rv,thetav

if(level > 0) then
  do i=1,mxp
    do j=1,myp
      do k=1,mzp
        thetav(k,i,j)=theta(k,i,j)*(1.0+1.61*rv(k,i,j))/(1.0+rtp(k,i,j))
      enddo
    enddo
  enddo
else
  do i=1,mxp
    do j=1,myp
      do k=1,mzp
        thetav(k,i,j)=theta(k,i,j)
      enddo
    enddo
  enddo
endif

return
end

!=======================================================================

subroutine exhtend(mzp,mxp,myp,ia,iz,ja,jz,pi0,pc,rtp,theta,  &
   rv,pt,thvtend,thvadv)

use rconstants
use micphys, only: level

implicit none


integer mzp,mxp,myp,ia,iz,ja,jz,i,j,k

real, dimension(mzp,mxp,myp) :: pi0,pc,rtp,theta,pt,thvtend,thvadv,rv

if(level.gt.0)then
  do j=ja,jz
    do i=ia,iz
      do k=2,mzp-1
        pt(k,i,j) = pt(k,i,j) + rocv * (pi0(k,i,j) + pc(k,i,j))  &
          / ( theta(k,i,j) * (1.0+1.61*rv(k,i,j)) / (1.0+rtp(k,i,j)) )   &
          * (  thvtend(k,i,j) + thvadv(k,i,j)  )
      enddo
    enddo
  enddo
else
  do j=ja,jz
    do i=ia,iz
      do k=2,mzp-1
        pt(k,i,j) = pt(k,i,j) + rocv * (pi0(k,i,j) + pc(k,i,j))  &
          / theta(k,i,j) * ( thvtend(k,i,j) + thvadv(k,i,j) )
      enddo
    enddo
  enddo
endif

return
end


!=======================================================================
subroutine exadvlf(m1,m2,m3,ia,iz,ja,jz,izu,jzv,rtgu,fmapui,rtgv,fmapvi   &
         ,f13t,f23t,rtgt,fmapt,dxt,dyt,uc,dn0u,vc,dn0v,dn0,wc,pt,pc)
     
use mem_grid

implicit none

integer i,j,k,m1,m2,m3,jm,im,ja,jz,ia,iz,jzv,izu

real :: c1z,c1x,c1y

real, dimension(m2,m3) :: rtgu,fmapui,rtgv,fmapvi,f13t,f23t,rtgt,fmapt  &
               ,dxt,dyt
       

real,dimension(m1,m2,m3) :: flxu,flxv,uc,dn0u,vc,dn0v,dn0,wc,flxw,pt,pc

! Compute momentum fluxes flxu, flxv, flxw

do j = 1,m3
   do i = 1,m2
      do k = 1,m1
         flxu(k,i,j) = uc(k,i,j) * dn0u(k,i,j) * rtgu(i,j)  &
            * fmapui(i,j)
         flxv(k,i,j) = vc(k,i,j) * dn0v(k,i,j) * rtgv(i,j)  &
            * fmapvi(i,j)
      enddo
   enddo
enddo

if(itopo == 0) then
   do j = 1,m3
      do i = 1,m2
         do k = 1,m1-1
            flxw(k,i,j) = wc(k,i,j)  &
               * .5 * (dn0(k,i,j) + dn0(k+1,i,j))
         enddo
      enddo
   enddo
else
   do j = 1,m3
      jm = max(j-1,1)
      do i = 1,m2
         im = max(i-1,1)
         do k = 1,m1-1
            flxw(k,i,j) = wc(k,i,j)  &
               * .5 * (dn0(k,i,j) + dn0(k+1,i,j))  &
               + hw4(k) * ((flxu(k,i,j) + flxu(k+1,i,j)  &
               + flxu(k,im,j) + flxu(k+1,im,j)) * f13t(i,j)  &
               + (flxv(k,i,j) + flxv(k+1,i,j)  &
               + flxv(k,i,jm) + flxv(k+1,i,jm)) * f23t(i,j))
         enddo
      enddo
   enddo
endif

! Compute advection contribution to U tendency

do j = ja,jz
   do i = ia,izu
      c1z = 0.5 / rtgt(i,j)
      c1x = c1z * fmapt(i,j) * dxt(i,j)
      c1y = c1z * fmapt(i,j) * dyt(i,j)

      do k = 2,m1-1
         pt(k,i,j) = pt(k,i,j) - c1x / dn0(k,i,j) * (  &
              flxu(k,i,j)  &
               * (pc(k,i,j) + pc(k,i+1,j))  &
            - flxu(k,i-1,j)  &
               * (pc(k,i,j) + pc(k,i-1,j))  &
            - (flxu(k,i,j) - flxu(k,i-1,j)) * 2.* pc(k,i,j) )
      enddo
   enddo
enddo

do j=ja,jzv
  do i=ia,iz
      do k=2,m1-1
        pt(k,i,j)=pt(k,i,j) - c1y /dn0(k,i,j) * ( &
          flxv(k,i,j)  &
           * (pc(k,i,j)+pc(k,i,j+jdim))  &
         -flxv(k,i,j-jdim)  &
           * (pc(k,i,j)+pc(k,i,j-jdim))  &
           -  (flxv(k,i,j)-flxv(k,i,j-jdim))*2.*pc(k,i,j) )
      enddo
   enddo
enddo

do j=ja,jz
  do i=ia,iz
      do k=2,m1-1
        pt(k,i,j)=pt(k,i,j) - c1z * dzt(k) /dn0(k,i,j) * ( &
          flxw(k,i,j)  &
           * (pc(k,i,j)+pc(k+1,i,j))  &
         -flxw(k-1,i,j)  &
           * (pc(k,i,j)+pc(k-1,i,j))  &
           -  (flxw(k,i,j)-flxw(k-1,i,j))*2.*pc(k,i,j) )
      enddo
   enddo
enddo

return
end

!=========================================================================

subroutine excondiv(m1,m2,m3,ia,iz,ja,jz,izu,jzv,uc,vc,wc,pc,pt  &
    ,dxt,dyt,rtgt,rtgu,rtgv,f13t,f23t,fmapt,fmapui,fmapvi)

use rconstants
use mem_grid
use micphys, only : level

implicit none

integer :: m1,m2,m3,i,j,k,ia,iz,ja,jz,izu,jzv,im,jm

real :: c1z,c1x,c1y
real, dimension(m1,m2,m3) :: flxu,flxv,flxw,uc,vc,wc,pt,pc
real, dimension(m2,m3) :: rtgu,fmapui,rtgv,fmapvi,rtgt,fmapt,dxt,dyt  &
         ,f13t,f23t
                             
! Compute divergence
!-----------
  ! Prep Fluxes
!-----------

do j=1,m3
  do i=1,m2
    do k=1,m1
      flxu(k,i,j)=uc(k,i,j)*rtgu(i,j)*fmapui(i,j)
      flxv(k,i,j)=vc(k,i,j)*rtgv(i,j)*fmapvi(i,j)
    enddo
  enddo
enddo

if(itopo == 0)then
  do j=1,m3
    do i=1,m2
      do k=1,m1-1
        flxw(k,i,j)=wc(k,i,j)
      enddo
    enddo
  enddo
else
  do j=1,m3
    jm=max(j-1,1)
    do i=1,m2
      im=max(i-1,1)
      do k=1,m1-1
        flxw(k,i,j)=wc(k,i,j) &
           + hw4(k) * ( (flxu(k,i,j)+flxu(k+1,i,j) &
           +flxu(k,im,j) + flxu(k+1,im,j)) * f13t(i,j) &
           + (flxv(k,i,j)+flxv(k+1,i,j) &
           +flxv(k,i,jm) + flxu(k+1,i,jm)) * f23t(i,j)) 
      enddo
    enddo
  enddo
endif

do j=ja,jz
  do i=ia,izu
    c1x=fmapt(i,j)*dxt(i,j)/rtgt(i,j)
    do k=2,m1-1
      pt(k,i,j)=pt(k,i,j)     & 
        - c1x * ( flxu(k,i,j)-flxu(k,i-1,j) )   &
        * pc(k,i,j) * rocv
    enddo
  enddo
enddo

do j=ja,jzv
  do i=ia,iz
    c1y=fmapt(i,j)*dyt(i,j)/rtgt(i,j)
    do k=2,m1-1
      pt(k,i,j)=pt(k,i,j)  &
        - c1y * (flxv(k,i,j)-flxv(k,i,j-jdim) )  &
        * pc(k,i,j) *rocv
    enddo
  enddo
enddo

do j=ja,jz
  do i=ia,iz
    c1z=1.0/rtgt(i,j)
    do k=2,m1-1
      pt(k,i,j)=pt(k,i,j)  &
        - c1z * dzm(k) * (flxw(k,i,j)-flxw(k-1,i,j) )   &
        * pc(k,i,j) * rocv
    enddo
  enddo
enddo

return
end

![ML
subroutine prep_vt3d(mxyzp,dtlto2,vt3da,vt3db,vt3dc,uc,vc,wc)
implicit none
integer                :: mxyzp, i
real                   :: dtlto2
real, dimension(mxyzp) :: vt3da, vt3db,vt3dc,uc,vc,wc

do i = 1,mxyzp
  vt3da(i) = (uc(i) + uc(i)) * dtlto2
  vt3db(i) = (vc(i) + vc(i)) * dtlto2
  vt3dc(i) = (wc(i) + wc(i)) * dtlto2
enddo

return
end subroutine prep_vt3d
!ML]
!-----------------------------------------------------------------------
!srf-
subroutine get_true_air_density(mzp,mxp,myp,ia,iz,ja,jz)
 USE mem_basic  ,  ONLY: basic_g
 USE mem_micro  ,  ONLY: micro_g
 USE rconstants ,  ONLY: cp,p00,cv,rgas,cpi
 USE micphys    ,  ONLY: level
 USE mem_stilt  ,  ONLY: stilt_g
 USE mem_grid   ,  ONLY: ngrid
 implicit none
 integer, intent (IN) :: mzp,mxp,myp,ia,iz,ja,jz
 real, parameter :: c1 = cv/rgas, c2 = p00/rgas !c2 = p00*(cpi**c1)/rgas
 real, dimension(mzp,mxp,myp) :: b
 real c3
 !real :: c3 = c2 * (cpi**c1)
 c3 = c2 * (cpi**c1)
 
 if( level == 0 ) then
     b(:,:,:) = 1.
 else
     ! b      = (1 + rtp)/(1+1.61*rv)
     b(:,:,:) = (1. + basic_g(ngrid)%rtp(:,:,:))/(1. + 1.61*basic_g(ngrid)%rv(:,:,:))
 endif 

!- true air density

 stilt_g(ngrid)%dnp(:,:,:) = ( c3/basic_g(ngrid)%theta(:,:,:)  ) * b(:,:,:) * &
 
                             ( basic_g(ngrid)%pi0(:,:,:) + basic_g(ngrid)%pp(:,:,:)) ** c1
 

!print*,'dnp=',stilt_g(ngrid)%dnp(2,20,20),basic_g(ngrid)%dn0(2,20,20)

!mxyzp = mxp * myp * mzp
!do i = 0,mxyzp-1
!  b=( 1. + a(irtp+i) ) / ( 1. + 1.61*a(irv+i) )
!  if(level==0) b=1.
!  a(idens+i)= (c2/a(itheta+i)) *b* ( a(ipi0+i)+a(ipp+i) )**c1
!enddo
end subroutine get_true_air_density


